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Abstract 

The ammonia dimer (NH3)2 has been investigated using high-level ab initio quantum chemistry 
methods and density functional theory (DFT). The structure and energetics of important isomers 
is obtained to unprecedented accuracy without resorting to experiment. The global minimum of 
eclipsed C s symmetry is characterized by a significantly bent hydrogen bond which deviates from 
linearity by as much as 20°. In addition, the so-called cyclic Cih structure, resulting from 
further bending which leads to two equivalent "hydrogen bonding contacts" , is extremely close in 
energy on an overall flat potential energy surface. It is demonstrated that none of the currently 
available (GGA, meta-GGA, and hybrid) density functionals satisfactorily describe the structure 
and relative energies of this nonlinear hydrogen bond. We present a novel density functional, 
HCTH/407+, which is designed to describe this sort of hydrogen bond quantitatively on the level 
of the dimer, contrary to e.g. the widely used BLYP functional. This improved generalised gradient 
approximation (GGA) functional is employed in Car-Parr inello ab initio molecular dynamics 
simulations of liquid ammonia to judge its performance in describing the associated liquid. Both 
the HCTH/407+ and BLYP functionals describe the properties of the liquid well as judged by 
analysis of radial distribution functions, hydrogen bonding structure and dynamics, translational 
diffusion, and orientational relaxation processes. It is demonstrated that the solvation shell of the 
ammonia molecule in the liquid phase is dominated by steric packing effects and not so much by 
directional hydrogen bonding interactions. In addition, the propensity of ammonia molecules to 
form bifurcated and multifurcated hydrogen bonds in the liquid phase is found to be negligibly 
small. 



* On sabbatical leave from: Department of Chemistry, Indian Institute of Technology, Kanpur, 
India 208016. 
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I. INTRODUCTION AND MOTIVATION 



Wavefunction-based quantum chemistry methods as well as density functional theory 
(DFT) are now well established electronic structure techniques within Theoretical 
Chemistry The former is traditionally used to investigate small finite systems, i.e. 
"gas phase problems", whereas the latter had its first successes in solid state physics, 
i.e., in the framework of "periodic lattice problems". There is also ample evidence that 
wavefunction-based quantum chemistry methods are able to predict properties of molecular 
systems up to virtually any desired accuracy as a result of systematically improving 
the correlation treatment together with the wavefunction representation. However, the 
numerical complexity, as e.g. measured by the computer time needed for such studies, 
explodes easily in praxi either before reaching a sufficient level of accuracy or before treating 
systems that are sufficiently large. DFT methods, on the other hand, are known to lead to 
astonishingly good predictions for a wide variety of systems at moderate cost - of course 
dramatic failures are also well documented in the literature. When describing molecular 
solids, associated liquids, or large biosystems, noncovalent interactions such as hydrogen 
bonds and van der Waals interactions might play an important role. The ideal approach 
would be to use post Hartree-Fock methods converged to the basis set limit, which turns 
out to be computationally prohibitive, especially when used in conjunction with molecular 
dynamics. While simple DFT methods are unable to correctly describe van der Waals 
interactions [3, [IB], hydrogen bonds are known to be well within the capabilities of DFT 
approaches once suitable gradient-corrected functionals are devised. A pioneering step in 
the early nineties was the investigation of hydrogen bonds in the water dimer and 
in small water clusters 8] based on such functionals, which are of the generalised gradient 
approximation (GGA) type. Soon thereafter it was demonstrated that also the structure and 
dynamics of liquid bulk water can be described faithfully [2J. Now there is ample evidence 
that various properties of water can be computed quite reliably using GG A-tvpe functionals, 
see e.g. Refs. E, S Q 11, Q Q 14,13,0 Q Q Q, Q and Refs. H Q, Q Q for 
dimer calculations. But how about other hydrogen bonded systems? 



Among hydrogen-bonded dimers the ammonia dimer is peculiar since it features a rather 
unusual hydrogen bond, see Refs. (25, 2(| for an early and a very recent overview. Instead of 
an essentially linear hydrogen bond, similar to the one of the water dimer having a HO- • -O 
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angle of only 5.5° 27], the hydrogen bond in the ammonia dimer is significantly bent as 
depicted in Fig. \Q or d. Furthermore, close in energy is the so-called cyclic structure where 
an even more pronounced bending of the hydrogen bond produces two equivalent "hydrogen 
bonding contacts" according to Fig. As an interlude we stress that such strongly bent 
hydrogen bonds play an important role in many biological systems, in particular if the 
bending allows for more than one interaction of a donor with acceptor molecules which 
leads to so-called bifurcated or multifurcated hydrogen bonds 0, 0|. m carbohydrates 
over 25% of the hydrogen bonds turn out to be multifurcated [28[ and in proteins |30j about 
40 and 90% of the hydrogen bonds in /3-sheet and a-helix structures, respectively, are of 
this type. In this context the ammonia dimer could serve as a prototype system in order to 
investigate nonlinear hydrogen bonding on the level of a minimal model, i.e. a free dimer in 
the gas phase. 

Guided by analogy or chemical intuition, one might guess that the ammonia dimer 
possesses a "classical" quasi-linear hydrogen bond similar to other dimers such as those 
of water or hydrogen fluoride. Early calculations, see e.g. Refs. HQ, Q and references 
cited therein, actually assumed from the outset a perfectly linear bond similar to the one 
shown in Figs. or b or predicted a quasi-linear hydrogen bond as depicted in Figs, or 



d, while others II 34 



favoured strongly bent up to cyclic dimers similar to Figs. ^ 



e. This controversy got fueled by microwave measurements [37[ providing evidence for a 
quasi-rigid cyclic structure such as the one in Fig. The issue of a quasi-linear hydrogen 
bond vs. a cyclic structure remained controversial (see e.g. Refs. 0, Q) up to early 



nineties where both experiment and theory started to converge. Thanks to soi 



experiments 
is now consensus 



e Dotn experiment and tneorys 
H £,0 00000 



listicated 
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49] there 



and extensive computations [4£ 
that the ammonia dimer is a fluxional molecular complex with an 
equilibrium structure that is characterized by a bent hydrogen bond such as the one shown 



0,00, 



in Fig. Combining far-infrared spectra with static and dynamic calculations 
the H-N- ■ -N hydrogen bond angle (between the N-H bond of the proton donor and the 
N- • -N axis) was determined to be 26° in the equilibrium structure with a N- • -N distance 
of ~ 3.354 A. Note, however, that this "experimental" structure of the ammonia dimer 



was deduced within the rigid monomer approximation 



systematic bias. In fact, a recent study 



48, 



which might introduce a 



of the water dimer has shown that monomer 



flexibility is indeed important in determining its vibration-rotation-tunneling spectrum and 
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thus "experimental" structures derived from it. It is expected that such monomer flexibility 
would also be important for the ammonia dimer. Furthermore, the global potential energy 
surface is found to be very flat BGU the^cCs— was estimated to be 
about only 7 cm -1 higher in energy, which is qualitatively consistent with the observation 
that the dimer is a very fluxional complex. Most interestingly, there is now experimental 
support [52! advocating a more cyclic equilibrium structure for (NH 3 ) 2 in superfluid Helium 
droplets at 0.4 K. 



53, 



541 ] on the structure 



More recent wavefunction-based quantum chemical calculations 
and energetics of the ammonia dimer support in some aspects the earlier conclusions 
Q, Q|- Second-order M0ller-Plesset perturbation calculations (MP2) with medium-sized 
basis sets (for geometry optimization followed by MP4 single-point energy calculations) [53] 
resulted however in a significantly smaller H-N- ■ -N hydrogen bond angle of 12° together 
with a reasonable low barrier of 7.6 cm -1 . Most recently, another (valence-only) MP2 
study j^, this time near the one-particle basis set limit, was reported. At the MP2 limit, 
the hydrogen bond angle of 23° was Cose, to the —enflefl ^ue EQ, tat th e ^ 
of 24 cm -1 seems too high and the dissociation energy D c (into separated and optimized 
monomer fragments) was estimated to be 13.5 ± 0..3 kJ/mol. This clearly hints to the 
importance of higher-order correlation effects and core-valence corrections, in particular 
since they account for 0.4 kJ/mol in the case oftflewate.dfln.Qtwhieh^exeeeo, 
the error bar given in Ref. j^J]). DFT calculations of the ammonia dimer are both rather 
scarce [^.l^l^l^ and somewhat inconclusive as to what level of accuracy DFT methods 
describe this intermolecular interaction. In the following we will clearly demonstrate that 
none of the 14 widely used GGA, meta-GGA, and hybrid functionals checked by us lead to 
a satisfactory description of both the structure and relative energies of the ammonia dimer 
(despite large basis sets being used and basis set superposition error being corrected for). 
Thus, even at the level of the ammonia dimer there is still room for improvement based on 
purely theoretical grounds ! 

Thus far we considered the behavior of ammonia in the limit of forming a dimer, e.g. in 
the very dilute gas phase. The question arises, if or how its peculiar behavior concerning 
hydrogen-bonding in the dimer is present in its condensed phases? Indeed, based on the 
apparent stability and rigidity of the cyclic dimer in the gas phase J^tJ it was concluded 
in quite general terms some time ago that NH3 might well be best described as a powerful 
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hydrogen-bond acceptor with little propensity to donate hydrogen bonds [25| ; no te that now 
it is however recognized that the ammonia dimer is non-ri gid and non-cyclic 26] . The x-ray 
structure of crystalline ammonia at low temperatures 0J^l| features a staggered, directed 
hydrogen bonded network, with the monomers being arranged like the dimer in Fig. ^i. 
Furthermore, based on neutron powder diffraction experiments, it was proposed that solid 
ammonia has an "unusual shared-hydrogen bond geometry" |62| and, in particular, the 
high-pressure phase ND3-/F contains bifurcated hydrogen bonds with bending angles as 
large as about 30°. A most recent theoretical study |63J], however, did not support the 
presence of bifurcated hydrogen bonds in phase IV of solid ammonia. In this latter study, it 
was concluded through calculations of electron density between intermolecular nitrogen and 
hydrogen atoms that the apparent bifurcated hydrogen bond geometry in solid ammonia 
IV is actually a single hydrogen bond perturbed by a neighboring interacting atom. The 
structure of liquid ammonia has been studied by x-ray 6J| and neutron 65] diffraction for a 
long time. Only the total structure factor and the total radial distribution function could be 
extracted in these early experiments and the presence of hydrogen bonds in liquid ammonia 
was inferred from a comparison of the experimental structure factor with the known structure 
of the solid. More recently, a set of neutron diffraction experiments was performed 6(| on 
liquid NH 3 , ND 3 , and an isomolar NH3/ND3 mixture. The powerful isotopic substitution 
technique, in conjunction with sophisticated data analysis, allowed extraction of all partial 
radial distribution functions at 213 and 273 K. Contrary to earlier conceptions (as presented 
in many standard textbooks on that subject) the spatial arrangement of nitrogen atoms 
showed that no extended hydrogen bonded network exists in liquid ammonia. Nevertheless, 
some degree of hydrogen bonding was inferred from the temperature dependence of the N-H 
and H-H radial distribution functions. However, the hydrogen bond interaction in liquid 
ammonia proved to be much weaker than that in water and no clear hydrogen bond peak 
was observed in either N-H or H-H correlations, unlike the case of water. 

Due to the importance of liquid ammonia and in particular its metallic solutions, also 
many theoretical attempts have been made in past two decades to understand its liquid 
phase. Most of these studies are based on classical Monte Carlo (MC) or molecular dynamics 
(MD) simulations using empirical site-site 

H H 0, H H Q 

or ab initio pair 



potentials [7^ which all lead to essentially linear hydrogen bonds in the limit of treating the 
ammonia dimer. A recent pioneering study 0] investigated the structure of liquid ammonia 
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0, i-e. 



through Car-Parrinello ab initio molecular dynamics [75|,|76[, i.e. without resorting to any 
pre-determined (pairwise) potentials but using the full many-body interactions as obtained 
by "on the fly" DFT (BLYP) instead. It was found that the results of the N-N, N-H and H-H 
radial distribution functions are in better agreement with experimental results than previous 
ones relying on model potentials. Some important differences, however, still remained. For 
example, the results of this study were found to overestimate the hydrogen bonded structure 

as revealed by the N-H correlations and predict a narrower first peak of the N-N distribution 

I I 

when compared with the corresponding experimental results. Also, a later such study [77j 
of liquid ammonia containing an ammonium or an amide solute reported somewhat different 
results for the N-N and N-H correlations between the bulk ammonia molecules compared to 
those of Ref. J^|. Besides, the nature of the hydrogen bonds in the liquid phase is not fully 
understood yet, in particular in view of the experimental situation. As discussed before, 
unlike water, the ammonia dimer has a bent hydrogen bond in its equilibrium geometry. If 
this bent character of the dimer hydrogen bond has the same significant effects on the liquid 
as in solid phases j3] is still an unresolved issue, in particular since all pair potentials by 
construction yield a quasi-linear hydrogen bond as do existing density functionals. 

The outline of this paper is as follows. In Section II, we describe the details of the 
computational methodologies that we employ for the calculations of the dimer and the liquid 
phase. In Section III, we present results for important dimer structures and their energies 
based on highly accurate ab initio methods and compare these data to those obtained from a 
wide variety of currently available density functional families. Section IV focuses in detail on 
the development of the novel HCTH/407+ functional, on its performance in various systems, 
and on the improvement achieved on the dimer level. In Section V, we compare the results 
of ab initio MD simulations of liquid ammonia using HCTH/407+ to those obtained from 
the widely used BLYP functional for reference purposes, and to experimental data wherever 
available. We summarise and conclude in Section VI. 
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II. METHODS AND TECHNICAL DETAILS 



A. The HCTH Density Functional Family 

In this contribution, we develop a new generalised gradient approximation (GGA) 
functional - which will be termed HCTH/407H — that only involves the density p(r) and its 
gradient Vp(r). The latter has been shown to be a necessary ingredient for the description 
of hydrogen bonds within DFT, see for instance Refs. 0, Q, 0, 22, 13, Q- To the 
contrary of with the so-called hybrid functionals that involve "Hartree-Fock exchange" 
(such as e.g. B3LYP0), 

this advantage is not impaired by a prohibitive computational 
overhead relative to the local (spin) density approximation if used in the framework of 
plane wave pseudopotential calculations, see e.g. the arguments presented in Sect. 3.3 of 



Ref. 



761 ] . Thus, we opted for reparametrising the Hamprecht- Cohen- Tozer-Handy (HCTH) 



form, which has been described in detail elsewhere 80] and was originally used by Becke 
for a hybrid functional jsjj. In general, HCTH is a post-local spin density approximation 
(post- LSDA) functional, meaning that it factorises the LSDA functional forms (Flsda) 

E xc = E i = EE c m FLSDA,~,{p a ,Pp)u g y {xl,xl)dr , (1) 

7=x,c [T , T ,c Q/ 3 7 g=o 

here, u g denotes the perturbation from the uniform electron gas if Co i7 is unity. The sum 
over the integers q going from zero to m implies the power series defined through it 7 , with 
9-y a and rj la as fixed coefficients: 

The variable x a is closely related to the reduced density gradient 

4 = £§£. (3) 

This functional form was previously employed for reparametrisations resulting in 

a highly accurate GGA functional, HCTH/407|s3]. Although the above form obeys none 
of the exact conditions for slowly varying densities 83], and violates part of the scaling 
relations Q], it does obey the most crucial scaling relations and the Lieb-Oxford bound 
0,0] within the physically important region 

When employing the form in Eqn. 1 up to fourth order in m, we obtain 15 linear 
coefficients (because of having exchange, spin-like and spin-unlike correlation), which are 



easily parametrised by minimising an error function f2, which is constructed out of the sum 
of three components 



n = £ w «(iC^-- E w" fl ) a + £«'i.G 



<9X 

+ E ^ / + ^ - vf- s ) a pfjdr . (4) 



The three sums represent the mean-square deviations from our reference data of energies, 
gradients and exchange-correlation potentials for each molecule, respectively, of the result 
of a Kohn-Sham density functional calculation (denoted by the superscript "K-S"). In the 
case of the total energies reference data are high-level quantum chemistry results that are 
denoted by "exact". In the second sum the exact gradients (at equilibrium geometry) 
are zero by definition and thus do not appear in the formula. In the final term, we fit 
to the exchange-correlation potentials as determined by the Zhao-Morrison Parr(ZMP)- 
method j^, which are shifted by a constant k because of the effects of the quantum- 
mechanical integer discontinuity. The ZMP method has been proven to be an important 
aspect of the fit J3]. All these contributions need to be balanced using weights w, which 
have been determined and reported previously The weights w consist of a product of 
multiple partial weights making contributions for each molecule in order to ensure a balanced 
functional. 

We also have to consider the molecular set for which the new functional was determined. 
In general, the new HCTH/407+ functional was determined by fitting it to the "407 set" 
of molecules as used for HCTH/407 0|. This set is similar to the G3 set of molecules j^jj, 
with added inorganic molecules, but a smaller proportion of large organic species which we 
considered to be over-represented. This original training set was supplemented with data 
from the ammonia dimer, including non-equilibrium structures, as described in detail in 
Sect. En 



B. Quantum Chemical Calculations 

For the accurate determination of some points of the potential energy surface (PES) 
of the ammonia dimer, we used the W2 method jj^]. This is basically an extrapolation 
towards the full CCSD(T) basis set limit including relativistic (but not Born-Oppenheimer) 
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corrections. W2 energies are calculated at the CCSD(T)/A'PVQZ geometry, this notation 
meaning an aug-cc-pVQZ basis set for nitrogen and a cc-pVQZ basis set for hydrogen 9^1 . 
The MOLPRO package Q] was used for these calculations. This method is one of the 
most accurate standard ab initio methods currently available, with an average error of less 
the, 0, M/mol for the GH M of »*cuU. U fn , it is -nore reiiable „d 

accurate than the Gl, G2 and G3 methods (although computationally much more expensive) 
HQ- 
Two types of higher-order contributions were considered. Firstly, higher-order T 3 effects 



were assessed by means of full CCSDT[98] calculations using the Aces 1 1 program system 



Second ly, th e importance of connect ed quadruple excitations was probed by me ans o f 
CCSD(TQ) PJhJ and BD(TQ) 



101 



calculations using AcesII and Gaussian98 

respectively. 

For the DFT calculations, we employed the Cadpac suite of pr ogram s 11041 . and assessed a 
variety of GGA and meta-GGA functionals (namely RLYP|1 Ofij. PBEpHfiJ, HCTH/120 



HCTH/407 



PW91 



HCTH[24jk and BP86[107j; we also tested but do not report results from 



108| and HCTH/93|8i|) and hybrids (B3LYP 3> B97-l |8Ql a nd r-HCTH hybri 



:idQ; 



109j and PBEO 110]). In case of the density 



not reported are the results obtained from B97-2 
functional calculations we used an A'PVTZ basis set, which is reasonably close to the DFT 
basis set limit. 

For the DFT values, which were obtained from using finite basis sets without 
extrapolati on a s done in the W1/W2 methods, we employed the counterpoise (CP) 
correction lllll to account for the basis set superposition error (BSSE). In the fitting 



procedure, the TZ2P basis set 



1121 was utilised in addition to the A'PVTZ basis set. 



C. Ab Initio Molecular Dynamics Simulations 



The ab initio molecular dynamics simulation s we re performed by means of the Car- 



Parrinello method 



75 



and the CPMD code 



76 



113j . A simple cubic box of 32 ammonia 



molecules with a box length of 11.229 A was periodically replicated in three dimensions 
and the elec tron ic structure of the extended system was represented by the Kohn-Sham 
formulation |ll4| of DFT within a plane wave bas is. The core electrons were treated via the 



atomic pseudopotentials of Goedecker et al. 



115] and the plane wave expansion of the KS 
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orbitals was truncated at a kinetic energy of 70 Ry. A fictitious mass of /x=875 a.u. was 
assigned to the electronic degrees of freedom and the coupled equations of motion describing 
the system dynamics was integrated by using a time step of 5 a.u. As usually done, see e.g. 
Refs. [3, 0, Q , the hydrogen atoms were given the mass of deuterium in order to allow for 
a larger molecular dynamics timestep which also reduces the influence of (the neglected) 
quantum effects on the dynamical properties. 

The ab initio MD simulations have been performed using the HCTH/407+ and 
BLYP |l05^ functionals. We used the BLYP functional in addition to the new functional 



because the former has been shown to p rovide a good description of hydrogen bonded liquids 
such as water 0, [ll] , methanol and also ammonia 0, 0| • Thus, it is worthwhile to 
compare the results of the two functionals for various structural and dynamical properties of 
liquid ammonia. The initial configuration of ammonia molecules was generated by carrying 
out a classical molecular dynamics simulation using the empirical multi-site interaction 
potential of Ref. jzjj. Then, for simulation with the HCTH/407+ (BLYP) functional, we 
equilibrated the system for 8.75 (10.1) ps at 273 K in NVT ensemble using Nose-Hoover 
chain method and, thereafter, we continued the run in NVE ensemble for another 9.30 
(9.34) ps for calculation of the various structural and dynamical quantities; the average 
temperatures of these microcanonical runs were about 275 (252) K. We note that the size of 
the simulation box used in the present simulations corresponds to the experimental density 
of liquid ammonia at 273 K which is 2.26 x 10~ 2 molecules/ A 3 



III. THE AMMONIA DIMER: STRUCTURE AND ENERGETICS 

In order to briefly validate the W2 method for some representative hydrogen bonded 
systems, (HF)2 and (H20)2, we compared the W2 results to MP2 basis set limit corrections 

, Ill7| , whereas the reference data was additionally 
empirically refined by scaling the calculations to certain quantum energy levels. Moreover, 
to get an estimate of the accuracy of the basis set extrapolation, we compared with the more 
cost-effective alternative of W2, which is Wl [92]. When calculated at the same geometry 
as W2, it usually yields a fair estimate as to how well the basis set limit is obtained in the 
extrapolation. In Table HJ the dissociation energies for the HF and H2O dimers are compared 
to best values from the literature, and displayed along with the results for the (NH3)(H20) 
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complex. All the Wl and W2 results are within the stated error margins of the reference 
data. While the dissociation energy of the hydrogen fluoride dimer is exactly the same, that 
of water dimer is slightly lower than predicted in reference In the latter case, relativistic 
corrections were not included in their best estimate. This will reduce the energy somewhat; 
however we would expect a geometry relaxation to lower the energy. The authors of Ref . |5j| 
used a CCSD(T)/aug-cc-pVTZ geometry and estimated the geometry relaxation, whereas 
the W2 calculation is carried out at an CCSD(T)/A'PVQZ geometry. Nevertheless, the 
results show only 0.03 kJ/mol difference between the Wl and W2 method for the HF dimer 
and 0.09 kJ/mol for the water dimer, indicating the high accuracy of these methods for such 
interactions. Finally, the (H20)(NH3) dissociation energy is pr edicted to be 26.8 kJ/mol by 
W2, compared to 27.4 kJ/mol determined by experiment |ll8j |: unfortunately no high-level 
ab initio data were found in the literature for this mixed dimer. 

Let us now focus on the ammonia dimer. Five different structures on the dimer PES can 
be considered important. The global minimum found on the CCSD(T)/A'PVQZ surface, as 
shown in Figured, has an HNN angle of 19.9° with eclipsed hydrogens and C s symmetry. 
Its staggered counterpart, displayed in Figure ^1, has an HNN angle of 13.1°. The two 
completely linear structures with the HNN angle being 0° (Figures ^ and ^d) are salient 
points and the C 2 h structure (Figured^) is a transition state. In Table HT) Wl and W2 results 
for the five structures are summarised, with column two and four (with BSSE) corresponding 
to the regular W2 method without the cp- correct ion. Because of the low barrier between the 
C s minimum and the C2h structure, we additionally calculated a counterpoise-corrected W2 
barrier. As to the accuracy of the electron correlation methods underlying the W2 method, 
we can make the following observations: 

• quasiperturbative connected triple excitations, i.e., the CCSD(T) - CCSD difference, 
only contribute 0.5 cm -1 to the C2/1 barrier with CP correction and 0.7 cm -1 without; 

• this has also been verified by CCSD(TQ) [BD(TQ)] calculations, where the 
contribution of the perturbative connected quadruples lowers the barrier by 0.2 cm -1 
[0.3 cm- 1 ] with the A'PVDZ basis set; 

• higher-order connected triples (i.e., CCSDT-CCSD(T)) additionally lower the barrier 
by 0.2 cm _1 when using A'PVDZ. This leads to the conclusion that not only has 
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basis set convergence been achieved, but also that higher-order excitations beyond 
CCSD(T) will contribute less than 1 cm -1 ; 

• further optimisation of the C s geometry at the MP2/ATV5Z level lowers the binding 
energy by only 0..3 cm" 1 (0.004 kJ/mol) relative to MP2/A'PV5Z//MP2/A'PVQZ; 

• we do not expect relaxation of the geometry at the core-valence level to make a 
significant contribution, si nce e ven for the monomer the core-valence correction to 
the geometry is very small |l lo| . 

Still, the difference between the BSSE-corrected and uncorrected results is rather large. 
Based on these results, we suggest that the C s to C 2 h energy difference is 3.5 ± 3 cm -1 , 
and we would expect the barriers between the linear and minimum geometries to be 
equally well described by W2. For the slightly improved CCSD(T) geometry, by correcting 
CCSD(T)/A'PVQZ - MP2/A'PVQZ + MP2/A'PV5Z, we predict the HNN angle to be 20.7 
degrees. As for the C 2 h barrier, the interaction energy is raised by 0.002 kJ/mol when 
going from CCSD(T) to CCSDT, showing that the quasiperturbative triples treatment (and 
the quadruples estimation as well) is quite accurate. Using an A'PVDZ basis set, the quasi- 
perturbative quadruples treatment lowers the binding energy by 0.092 kJ/mol for the coupled 
cluster and 0.087 kJ/mol when using the Brueckner-Doubles method. With the geometry 
relaxation contributing less than 0.005 kJ/mol at MP2 level and the estimate of the higher 
order excitations, we predict the dissociation energy to be 13.1 ± 0.2 kJ/mol, since the main 
error is likely to come from the quadruple excitations. 

The results obtained using a selection of currently published functionals are shown in 
Table HH For comparison, the CCSD(T), CCSD, MP2 and HF values are given. It should 
be noted, however, that BSSE for the coupled-cluster values is more than twice as high as 
for DFT (0.9 compared to 0.4 kJ/mol), and therefore we would expect the DFT calculations 
to be closer to their respective basis set limits than the ab initio methods. Nevertheless, 
their BSSE is still larger than that at the HF level. As expected, the CCSD(T) and W2 
numbers are reasonably close to each other. Most density functional methods reproduce the 
dissociation energy as listed in the first column, reasonably well. In general, the GGA and 
meta-GGA functionals tested give a good description of this interaction underestimate and 
PW91 overestimating the dissociation energy. The tested hybrid functionals are even more 
accurate, yet all of them underestimating the dissociation energy. 

13 



However, if we consider the energy difference between the staggered and eclipsed 
conformation of the NH 3 -dimer (see the data in the second column of Table ITTTj) . it becomes 
obvious that available functionals are incapable of reproducing the effect of a bent hydrogen 
bond. It is found that density functionals typically underestimate the difference between 
the local and global minima by at least a factor of two. We note here that functionals 
other than those listed in our tables (such as e.g. PW91, B97-2, PBEO) yield similar results 
which are therefore not displayed. On the other hand, all other wavefunction-based ab initio 
methods considered are able to reproduce this effect. This becomes even more apparent 
when considering the C2h transition state. Here, the coupled-cluster differences are close to 
3 cm -1 , whereas HF renders 24 cm -1 and the DFT values are between 60 and 95 cm -1 with 
the exception of HCTH/407 (27 cm" 1 ) and the r-HCTH functionals (100 and 140 cm- 1 ). In 
turn, the other barrier is underestimated by 30-60% and both linear structures are stabilised 
relative to the global minimum. The same trends are seen in Table HVl where the HNN angles 
of the two minimum structures are shown. Here, available density functionals underestimate 
the bending of the hydrogen bond by almost a factor of two for the global minimum and by 
about 30% for the local (staggered) minimum, although the hydrogen bond distance itself 
is reasonably well reproduced. Both r-HCTH functionals, while they perform well with 
their large 'training' fit set, are the worst performers here. On the other hand, HCTH/407 
is the only functional which is somewhat capable of reproducing the effect, outperforming 
Hartree-Fock. 

In summary, although standard density functionals do yield a correct dissociation energy 
and a slightly nonlinear hydrogen bond, they completely overestimate the barrier towards 
the symmetrical C^ structure, thereby preferring a much more linear hydrogen-bonding 
configuration than our accurate reference methods. Although the energy differences might 
be viewed as intrinsically too small to be accurately rendered by DFT, the fact that even HF 
gives better results is somewhat disturbing. Nevertheless, the observation that Hartree-Fock 
and HCTH/407 are somewhat capable of properly predicting the structure and energetics 
of the ammonia dimer system suggests that the problem at hand is not beyond the grasp of 
density functional methods as such. 
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IV. THE HCTH/407+ FUNCTIONAL: CONSTRUCTION AND PERFOR- 
MANCE 



Considering that standard density functionals were shown to be unable to describe the 
energy difference between linear and bent hydrogen bonds and do not capture the large 
bending angle of the hydrogen bond, the most obvious remedy would be to change the 
functional form itself. However, if all tested functionals give such a bad description, this 
is unlikely to solve the problem. The culprit in this case, therefore, has to be the set of 
molecules which the HCTH functionals are fit to, insofar as only linearly hydrogen-bonded 
systems have been included. Nevertheless, including the minimum structure and energetics 
of the ammonia dimer determined in the previous section into the set probably only partly 
solves this problem. The main problem, as established before, is the energy difference 
between the linear and bent hydrogen bonds. In addition, the TZ2P basis set commonly 
used in the fit set might not be appropriate. Even when employing the A'PVTZ basis 
set mentioned above, the basis set superposition error for the dimerisation energy is 29 
cm -1 , contributing 8 cm" 1 to the energy difference between the minimum and symmetric 
structure when using the B97-1 functional. In comparison, using the TZ2P basis set, the 
BSSE is 92 cm -1 , contributing 33 cm -1 to the energy difference. Clearly, the TZ2P basis set 
is inappropriate for this problem. The above numbers indicate that a counterpoise procedure 
has to be employed within the fit in order to get a good estimate. Note that even when 
using an extrapolation to the full basis set limit (in the case of W2), the difference between 
the counterpoise corrected and uncorrected values is already 3 cm -1 . 

All of the above issues were considered in the fit.. Thus, we fit to all points shown in 
Figure employing an A'PVTZ basis set with counterpoise correction in the fit. This is 
the first time we are fitting to non-equilibrium structures, and we only fit to their energy 
differences and exchange-correlation potentials of these. Due to the energy difference of only 
a couple of wavenumbers, an additional weight has to be employed for any change in the fit 
to be noticeable. Thus, equation (3) now reads 




•K-S\ 2 



<K-S 



jpBSSE\2 , 



m 



m 



l,X 




(5) 



15 



with 



rpBSSE rpghost, donor , rpghost, acceptor rpmonomer, donor rpmonomer, acceptor /r?\ 

m m ' m m m ' V u / 

Since we do not want to fit to a given value of basis set superposition error, it is included in 
E^~ s . The weights now also include the additional weights for the potential energy 

surface with the other weights defined as before 0|. Thus, w m is composed of 

l PES = 750 x weight con/idence x weight aiom x weight dimer x weight PES . (7) 



w 



When determining the weights for the new members of the fit set, we need to consider 
both their A'PVTZ optimised energies as well as their A'PVTZ single-point energies at 
the CCSD(T)/A'PVQZ optimised geometries. The latter is necessary since we can only fit 
to single points. As a result of the small energy differences, the frozen-geometry single- 
point energies and optimised energies differ by a significant amount. The development 
procedure is similar to the one employed for the HCTH/407 functional, and the set includes 
numerous atomisation energies, proton affinities, ionisation potentials, electron affinities, 
and dissociation energies of hydrogen-bonded dimers and transition metal complexes 82]. 
First, we fit to a subset of the 407 molecules in order to get close to the global minimum, then 
we fit to the full set. The latter subset (small set) is the G2-1 set plus two hydrogen-bonded 
systems (H 2 0)2 and (HF) 2 , and the additional structures on the potential energy surface, 
making 153 systems in all (the 147 of the HCTH/147 functional plus the five PES points, in 
addition to the NH 3 dissociation energy at the TZ2P basis). For each point on the PES, the 
counterpoise correction was employed. The results using the newly determined functionals 
as a function of the weights have been determined in both the first fit (to 153 systems) and 
the second fit (to 413 systems). Thus, we obtained new functionals with weights ranging 
from 0/0 (which corresponds to the HCTH/407 functional) to 200/200. The first weight is 
the one used in the subset, while the second corresponds to the large fit. The importance 
of the weights in the first fit can be seen from the differences by using weight pairs of 
30/80 and 80/80 - they differ by about a factor of two in the energy difference between 
the staggered and eclipsed geometry. Hence, the minimum obtained by the second fit will 
of necessity be close to the global minimum obtained by the first one. When increasing 
these weights, we approach the reference values as expected, with the 200/200 functional 
rendering the closest values to the reference method. While not explicitly included in the fit, 
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the BSSE also increased from 0.444 kJ/mol to 0.866 kJ/mol. This resembles the tendency 
of the functional to have a larger intermolecular distance than the reference methods when 
increasing the weights, hence entering a regime where the BSSE increases more rapidly than 
for the optimised structures. In comparison, the largest BSSE observed by the functionals 
tested in Table IIHI is 0.48 kJ/mol. However, for all the calculations in the fit we have 
performed only DFT single-point energy calculations at the CCSD(T) optimised geometries. 

Another, more appropriate, point of comparison consists of the optimised structures of 
all the newly obtained functionals, as done in Table IIHI for the contemporary functionals. 
Since the energy differences in general are very small, these results differ significantly from 
the ones obtained by the single-point energies calculated at the coupled-cluster geometries. 
In addition, the C2h structure becomes the global minimum for the functionals with a weight 
of 80 or larger in the first fit. Based on these results, the 30/80 functional appears to be 
the most reliable for this interaction, and hence will be denoted as the new HCTH/407+ 
functional. This new functional (compare to Table IIII)) now renders a dissociation energy 
of 13.18 kJ/mol for the ammonia dimer, in much closer agreement to the reference values 
than any other method employed. The difference between the C s minimum and the C2/1 
structure is 4.0 cm -1 for HCTH/407+, compared to our best estimate of 3.5 cm" 1 in Table 
ITT1 The staggered C s structure lies 18 cm -1 above its eclipsed counterpart, which is about 
five wavenumbers higher than our best estimate but still within its error bars. The difference 
to the linear eclipsed and linear staggered structures are slightly underestimated at 51.5 and 
53.5 cm -1 respectively, but this is still a vast improvement over the other functionals. In 
Figure El we compare the energies of the new HCTH/407+ functional to W2 and one of the 
most commonly used GGA functionals (BLYP). The BSSE for the relaxed C s structures does 
not increase as rapidly as for their CCSD(T) optimised counterparts. It increases from 0.43 
(HCTH/407) to 0.45 (HCTH/407+) to 0.61 kJ/mol (weight of 200), and does not yield as 
unreasonably large values as the single-point calculations at the CCSD(T) geometries. The 
hydrogen bond lengths of the C s and C2/1 structures increase with larger weights, emphasising 
the importance of this analysis. Generally, all HCTH functionals are found to slightly 
overestimate this hydrogen bond length. The HNN angle of the HCTH/407+ functional is 
17 for the minimum geometry. This is within the accuracy that can be expected from 
contemporary density functionals (± 2°), showing a significant improvement over all other 
density functionals in Table HVl The results for all functionals obtained with different weights 
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(from whi ch w e determined the HCTH/407+ functional) are given in the supplementary 



material. 



120] 



The parameters which define the HCTH/407+ functional, compared to HCTH/407, are 
given in Table IVl As might be expected, they differ only slightly from those for the standard 
HCTH/407 functional; only the correlation and higher-order coefficients have significantly 
been affected by the change.. Nevertheless, the new functional should now be able to describe 
non-directed hydrogen bonds better than HCTH/407. On the level of the ammonia dimer, 
the improvement with respect to other density functionals becomes clearly visible in Fig. |21 
see Table ITTTI for the corresponding numbers. It is seen that HCTH/407+ describes the 
relative stability of all isomers much better with respect to the W2 reference data than BLYP, 
which was selected as a prominent representative of the GGA family. Most importantly, the 
dramatic failure of BLYP to capture the stability of the cyclic C2/1 isomer relative to all 
other isomers is cured. According to HCTH/407+, the two C s structures as well as the 
C2h structure are essentially degenerate, which is in accord with the reference data.. Apart 
from the energetics, also the structure of the ground state of the ammonia dimer, i.e. the 
eclipsed C s structure, is dramatically improved, see Table H^l for details. In particular the 
hydrogen bond angle HNN is now much larger and thus closer to the reference value than 
any other density functional method; note that also the functionals PW91, B97-2, and PBEO 
were considered. The same is true for the HNN angle in the staggered variant of the C s 
symmetric structure. However, it is also clear from this table that the N- • -H distance is 
clearly overestimated: HCTH/407+ yields a distance of about 2.5 A instead of around 
2.3 A as required. This seems to be typical for the HCTH family as both HCTH/120 and 
HCTH/407 yield similarly large values exceeding 2.4 A (and HCTH/93 leads to 2.614 A). 
In addition, the same trend of producing too long hydrogen bonds is found for the C^h 
structure. This deficiency is corrected by both r-HCTH and r-HCTH hybrid functionals. 
Here it should be noted that both BLYP and PBE describe the hydrogen bond length quite 
well, however at the expense of making it much too linear and grossly disfavoring the cyclic 
structure energetically. 

Furthermore, this significant improvement in describing non-directed hydrogen bonds, 
however, comes at the expense of a slightly increased error for the other molecules in the 
fit set. Table EU shows the errors of the HCTH/407 and HCTH/407+ functionals to the 
fit set, summarising these as RMS energy and geometry errors, together with the errors for 
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the energies, geometry shifts and frequency shifts of nine hydrogen-bonded systems. The 
values for B3LYP, BLYP and BP86 are included for comparison. The hydrogen-bonded 
systems are the (HF) 2 , (HC1) 2 , (H 2 0) 2 , (CO)(HF), (OC)(HF), (FH)(NH 3 ), (C1H)(NH 3 ), 
(H 2 0)(NHs) and (H30 + )(H 2 0) complexes, but not the NH3 dimer itself. The performance 
of other functionals over all the sets has been discussed elsewhere 0, 1^, H |. Overall, 
errors for the HCTH/407+ functional are quite similar to HCTH/407, which in turn means 
that they are still excellent compared to first-generation hybrid functionals like B3LYP 
and especially GGA functionals like BLYP or BP86. Although we expect the geometries 
to be still slightly worse than B3LYP (although the gradient error is comparable), the 
energetic properties are better described, as are the hydrogen bonds. Whereas the errors 
of the HCTH/407+ functional for the fit set have barely changed compared to HCTH/407, 
the accuracy of the dissociation energies of the nine directed hydrogen bonds has actually 
decreased by about 20%. Nonetheless, the accuracy of the frequency shifts of the H-X bond 
stretches has increased, indicating that we seem to have improved the description of the 
potential energy surface of these complexes. 



V. LIQUID AMMONIA: STRUCTURE, DYNAMICS, AND HYDROGEN BOND- 
ING 



A. Radial Distribution Functions 



We have calculated the nitrogen-nitrogen, hydrogen-hydrogen and nitrogen-hydrogen 
radial distribution functions (RDF) from the atomic trajectories generated by ab initio 
molecular dynamics simulations; note again that we used the mass of the deuteron instead 
of that of the proton but we still use the symbol H for convenience. The results for both 
HCTH/407+ and BLYP functionals are shown in Fig. EJ In this Figure, we have also 
included the experimental results 0, Il2l| of the RDFs of liquid ammonia at 273 K; note 
that the underlying experimental technique is based on the assumption that the structure 
of all isotopically substituted systems is identical. For the nitrogen-nitrogen RDF, it is seen 
that at the short distances the BLYP gNN{f) is in better agreement with experimental 
results. The HCTH/407+ peak is located at a distance of about 0.15 A larger than 
the experimental peak position whereas the BLYP peak appears at a somewhat shorter 
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distance. We note in this context that the HCTH/407+ functional gives an overestimated 
equilibrium N-N distance for the isolated dimer compared to the reference CCSD(T) result 
as can be inferred from Table IIVI The larger N-N distance seems to be translated to the 
liquid phase configurations and we observe a somewhat larger most probable N-N distance 
for the HCTH/407+ functional than what is found in the experiments. Still, the overall 
shape of the first peak is represented reasonably well by both the functionals. When 
the RDFs are integrated up to their first minimum to obtain the coordination number 
in the first solvation shell, we obtained the values of 13.2 and 12.1 for HCTH/407+ and 
BLYP functionals compared to an experimental value of 12.75. Clearly, the agreement with 
experimental results is found to be rather good. 

Considering the results of the H-H and N-H RDFs, the positions of the calculated 
intramolecular H-H and N-H peaks for both the functionals agree rather well with 
experimental results as expected. For example, the liquid phase N-H bond length and 
intramolecular H-H distance for the HCTH/407+ (BLYP) functionals are 1.03 (1.03) A and 
1.63 (1.65) A which can be compared with the experimental values of 1.04 A and 1.60 A, 
respectively. Also, these intramolecular distances in the liquid phase are very close to their 



gas phase values. This was also noted earlier for the BLYP functional |59|. Because of the 
quantum disperson effects which are present in real liquid but not in the current simulated 
systems, the experimental intramolecular peaks are found to be broader than the theoretical 
peaks. The somewhat higher value of the experimental H-H correlation at the first minimum 
at around 2 A can also be attributed to such quantum dispersion effects. 

The intermolecular H-H peaks are found to be well represented by both functionals. The 
intermolecular part of the N-H RDF shows the presence of both H-bonded and non H-bonded 
H atoms in the first solvation shell. The shoulder up to about 2.7 A can be attributed to 
hydrogen bonded H atoms, while the more pronounced peak at the larger distance at around 
3.75 A corresponds to the hydrogen atoms in the solvation shell which are not hydrogen 
bonded. Again, both functionals are found to describe well the non-hydrogen bonded first 
solvation shell peak when compared with the corresponding experimental result. For the 
hydrogen bonded part, however, the HCTH/407+ result is found to exhibit a less pronounced 
shoulder than what is found in the experiments and for the BLYP functional. This is 
probably a result of the slight overestimation of the hydrogen bond length by HCTH/407+ 
already on the level of the dimer, see Table IIVI Considering both hydrogen bonded and 
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non-hydrogen bonded H atoms, a sphere of radius 5.3 (5.. 2) A around a central N atom is 
found to contain 34.1 (34) intermolecular H atoms for the HCTH/407+ (BLYP) functional. 
The radius of the sphere is set to 5.3 (5.2) A because the first minimum of the intermolecular 
N-H RDF is located at this distance. The corresponding radius for the experimental N-H 
RDF is 5.0 A, yielding the experimental value of 34.8 for the above coordination number. 
A more detailed analysis of the distribution of hydrogen bonds in liquid ammonia is given 
in the following subsection, and the dynamical aspects of hydrogen bonds are considered in 
Sect. IVDl 



B. Distribution of Hydrogen Bonds 

The analysis of the hydrogen bond (HB) distribution among ammonia molecules is based 
on a calculation of the fractions P n of ammonia molecules that engage in n hydrogen 
bonds and the average number of hydrogen bonds that an ammonia molecule donates and 
accepts [2 



liquids 



Follo wing previous work on liquid amm onia, water and other hydrogen bonding 
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71, 72, 122, 123 



124, 



125 




Il28l | , we have adopted a geometric definition 



of the hydrogen bonds where two ammonia molecules are assumed to be hydrogen bonded 
if they satisfy the following configurational criteria with respect to N-N and N-H distances 

R (NN) < R (NN) ) 

R (NH) < R (NH) (g) 

and the values of these distance cut-offs are usually determined from the positions of 
the first minimum of the corresponding (intermolecular) RDFs. We have used a similar 
procedure where the N-N RDF gives a value of R™ =5.25 (5.1) A for the N-N cut- 
off for HCTH/407+ (BLYP) functionals. However, the procedure can not be applied 
unambiguously to determine the N-H cut-off because no clear minimum that separates the 
hydrogen bonded and nonbonded H atoms is found in the N-H RDF. In the latter case, 
the shoulder which is assigned to hydrogen bonded H atoms merges into the broad and 
more intense peak corresponding to the H atoms of the first solvation shell which are not 
hydrogen bonded. Since this shoulder - which corresponds to the hydrogen bonded H atoms 
- extends up to about 2.7 A in the experimental and BLYP RDFs we have taken this 
distance of 2.7 A as the cutoff N-H distance. For HCTH/407+ also, we have used the same 
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cut-off of the N-H distance although the hydrogen bonded shoulder is less pronounced for this 
functional. Still, variation of R^ H within reasonable bounds does not change qualitatively 
the result of this analysis. 

In Fig. HJ we have shown the distribution of donor and acceptor hydrogen bonds for both 
HCTH/407+ and BLYP functionals. In this figure, panels (a) and (b) show the fraction 
of ammonia molecules that accept or donate n hydrogen bonds and panel (c) shows the 
fraction of hydrogen atoms that donate n hydrogen bonds. The majority of the ammonia 
molecules is found to accept one hydrogen bond and to donate one. Also, the fraction of 
hydrogen atoms donating two or more hydrogen bonds is found to be negligibly small, which 
means that the bifurcated or multifurcated hydrogen bonds are practically absent in liquid 
ammonia. The average number of donor hydrogen bonds per ammonia molecule, which is 
also equal to the average number of acceptor hydrogen bonds per molecule, is found to be 
0.93 (1.34) for the HCTH/407+ (BLYP) functional. Given that there are as many as about 
13 ammonia molecules in the first solvation shell, we conclude that non-hydrogen-bonded 
interactions, i.e. packing or steric effects, are crucial in determining the solvation behavior 
of ammonia in the liquid state. This is in sharp contrast to the case of liquid water where 
the solvation structure of a water molecule is determined primarily by hydrogen bonding 



interactions 



129|. 



The ammonia dimer has a strongly bent hydrogen bond in its equilibrium geometry and 
it would be interesting to investigate whether the hydrogen bonds remain bent in the liquid 
phase. We have carried out such an analysis by calculating the distribution function of the 
cosine of the hydrogen bond angle 8, which is defined as the N- ■ • N-H angle of a hydrogen 
bonded pair and the results are shown in Fig..|2k for both the functionals. It is seen that, for 
both the functionals, the probability distribution P(cos8) is peaked at cos# = 1 indi cating 
that the preferred hydrogen bond geometry is linear in the liquid phase, see Ref. jl3o| . The 
peak is somewhat more sharp and narrower for the BLYP functional, which means a stronger 
preference for the linear hydrogen bonds for this functional as compared to that for the 
HCTH/407+ functional. This preferential presence of linear hydrogen bonds in the liquid 
phase is indeed an interesting result given that the dimer has a bent hydrogen bond in the 
gas phase. We have also investigated the distribution of N- • • N-H angle for H atoms that 
belong to nearest neighbours but are not hydrogen bonded, i.e. the N- ■ -H distance is greater 
than 2.7 A. The results of these angular distributions are showninFig.Ebfor2.7A < R^) 
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< 3.7 A and in Fig. Efc for 3.7 A < R^) < 4.7 A. A wide range of values are found for the 
N- • • N-H angle with very little or no preference for a particular orientation. These findings 
are indeed consistent with the above analysis that the first solvation shell in liquid ammonia 
is dominated by simple packing rather than by directional hydrogen bonding. 



C. Self Diffusion and Orientational Relaxation 



The translational self diffusion coefficient of an ammonia molecule in the liquid state is 
calculated from the long-time limit of the mean-square displacement (MSD) 

= <[r( f )-r(0)]'> 

where r(t) is the position of the center of mass of a molecule at time t and the average 
is carried out over the time origin for autocorrelation and over all the molecules as usual. 
The results of the diffusion coefficients, which are obtained from a least square linear fit 
of the simulation data excluding the initial ballistic regime up to 0.5 ps, are included in 
Table (VII). This Table also contains the experimental values of the diffusion coefficie nts o f 
deuterated liquid ammonia ND3 at 275 and 252 K that are found from Eq. (2) of Ref . |l3l| , 
which was shown to fit the experimental data over a rather wide temperature range very 
well. The agreement of the calculated diffusion coefficients with the experimental result is 
found to be reasonably good for both functionals. In particular it is gratifying to see that 
the experimental value is underestimated because quantum effects on nuclear motion, which 
are neglected in the present study, tend to increase the diffusion in a liquid. 

The self-orientational motion of a mmo nia molecules is analyzed by calculating the 
orientational time correlation function |l32j |. Cf (t), defined by 

ca(t) _ < Pi(e a (t) • e"(0)) > 

where Pi is the Legendre polynomial of order / and e a is the unit vector which points along 
the a-axis in the molecular frame. Here, we have studied the time dependence of Cf (t) 
for I = 1,2 and for three different e a , the unit vectors along the molecular dipole axis, 
an N-H bond and an intramolecular H-H axis. The geometric dipole vector of an ammonia 
molecule is calculated by assigning partial charges to N and H atoms as given by the classical 
model of Ref.jzj]. For the orientational relaxation of unit vectors along the N-H and H-H 
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axes, the results are averaged over three such intramolecular axes of each type. After an 
initial transient non-exponential decay, the relaxation becomes diffusional and C^(t) decays 
exponentially. The orientational correlation time, r™, is defined as the time integral of the 
orientational correlation function 







dt C?(t) . (11) 



We note that the orientational relaxation of molecular vectors containing H atoms are usually 
measured by NMR relaxation experiments, which yield the Fourier transform of a correlation 
function for I = 2 such that the integrated correlation time r-f is measured. 

In Fig. El we depict the results of the orientational correlation function C"(t) for I = 1,2 
and a=dipole, N-H and H-H. The corresponding results of the orientational correlation times 
are included in Table (VII). We have also included the experi menta l result of T2 of deuterated 



ammonia ND% as given by the fitted Arrhenius law of Ref. |1 331 ] . Since this experimental 
relaxation time was obtained by NMR relaxation experiments, it is likely to correspond to 
the relaxation of N-H/H-H vectors rather than the dipole vector. The BLYP functional is 
found to predict a relatively slower orientational relaxation. Considering the quite small 
system size and short simulation trajectories, both functionals are found to describe the 
single-particle rotational dynamics of liquid ammonia reasonably well. 



D. Dynamics of Hydrogen Bonds 



The key dynamical quantity in the context of hydrogen bond dynamics is the mean 
hydrogen bond life time, thb- The fast librational motion of ammonia molecules makes 
apparent breaking and reformation of a hydrogen bond a very fast process and, depending 
on how these fast transient events are taken into account, one can obtain different values 
for the hydrogen bond (HB) lifetime. The rotational motion of ammonia molecules is also a 
rather fast process and can also contribute to the short-time dynamics of hydrogen bonds. 



For a p articular chose n definition of a breaking event, one can 



method 



128, 



134L Il3.j | or a time correlation function method 
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adopt a direct counting 
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1221, 136] 



to calculate the lifetime. In the present work, we have adopted the time correlation function 
approach which allows us to calculate the hydrogen bond lifetime in both cases: (a) The 
breaking of hydrogen bond occurs due to fast librational motion even though it may be a 
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transient event and not a 'true' breaking event and (b) these fast librational breaking events 
are ignored and only rotational and translational diffusional breaking on a longer time scale 
is considered as the 'true' breaking of a hydrogen bond. Also, in the present context, it is 
implied that 'fast librational breaking' also includes the effects of fast rotational motion that 
might contribute to the short time dynamics of the hydrogen bonds. 

In order to obtain the hydrogen bond lifetime for the two scenarios mentioned above, we 
calculate three time correl ation function s: one continuous and two intermittent hydrogen 
bond correlation functions 123 . 124 . 13(| . Before we define these time correlations, we first 
define two hydrogen bond population variables h(t) and H(t): h(t) is unity when a particular 
tagged pair of ammonia molecules is hydrogen bonded at time t, according to the adopted 
definition from Sect. IV 51 and zero otherwise. The function H(t) is unity if the tagged pair 
of ammonia molecules remains continuously hydrogen bonded from t — to time t given 
that the bond was formed for the last time at t = 0, and it is zero otherwise. We define the 
continuous hydrogen bond time correlation function SnBit) as 

S HB (t) =< h(0)H(t) > I < h > , (12) 

where < • • • > denotes an average over all hydrogen bonds that are created at t = 0. Clearly, 
Shb (t) describes the probability that a hydrogen bonded ammonia pair, which was created 
at t — 0, remains continuously bonded upto time t. It differs from the continuous hydrogen 
bond correlation function of Rapaport [136] in that we apply the condition of bond formation 
at time t = 0. The mean hydrogen bond lifetime the can then be calculated from the time 
integral 

Thb = / S HB (t)dt . (13) 
Jo 

In Fig. Et) we have shown the decay of SnB{t) for both functionals and the corresponding 
results of thb are included in Table VII. The hydrogen bond lifetime is found to be about 
0.1 ps for both functionals. We note that both fast librational and slower diffusional motion 
can contribute to the decay of Shb(P)- Since the librational motion occurs on a faster 
time scale and since the correlation function SnB{t) does not allow any reformation event, 
the dynamics of SnB{t) primarily reveal the dynamics of hydrogen bond breaking due to 
fast librational motion and hence the lifetime that is obtained from Eq.(12) corresponds to 
the average time over which a hydrogen bond survives before it 'breaks' due to librations. 
Again, this short time breaking may also include contributions from fast rotational motion 
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of ammonia molecules which also occurs at subpicosecond time scale. 

A different way to analyze the h ydrogen b ond dynamics is to calculate the intermittent 



hydrogen bond correlation function 



123, 



C HB {t) =< h{0)h{t) > I < h > , (14) 

where the average is now over all hydrogen bonds that were present at t — 0. Note that the 
correlation function CuB{t) does not depend on the continuous presence of a hydrogen bond. 
It describes the probability that a hydrogen bond is intact at time t, given it was intact at 
time zero, independent of possible breaking in the interim time. Clearly, bonds which were 
briefly 'broken' by fast librational motions would continue to contribute to the correlation 
function at later times and this leads to a much slower decay of CnB{t) at longer times. The 
results of CnB{t) are shown in Fig. 03 which clearly shows a two-phase relaxation of this 
correlation function. The initial fast relaxation corresponds to rapid breaking of hydrogen 
bonds due to librational motion and the slower relaxation after this initial transient decay 
corresponds to the breaking and reformation of hydrogen bonds due to both rotational and 
translational diffusion of ammonia molecules. 

After a hydrogen bond is broken, the two ammonia molecules can remain in the vicinity 
of each other for some time before either the bond is reformed or the molecules diffuse away 
from each other. We define NnB{t) as the time dependent probability that a hydrogen bond 
is broken at time zero but the two molecules remain in the vicinity of each othe r i.e. as 



nearest neighbors but not hydrogen bonded at time t.. Following previous work 123], we can 
then write a simple rate equation for the 'reactive flux' —dCnB/dt in terms of C##(t) and 
N HB (t) 

- = kC HB (t) - k'N HB {t) , (15) 

where k and k! are the forward and backward rate constants for hydrogen bond breaking. 
The inverse of k can be interpreted as the average lifetime of a hydrogen bond. 

The derivative of the intermittent hydrogen bond correlation or the 'reactive flux' of 
Eq.(14) is computed from the simulation results of Chb^) that are presented in Fig. [7)3. 
The probability function Nhb^) is also calc ulate d from the simulation trajectories through 
the following correlation function approach |l24| 

N HB (t) =< h(0)[l - h(t)]h'(t) > I < h > , (16) 
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where h'(t) is unity if the N-N distance of the pair of ammonia molecules is less than 
at time t and it is zero otherwise. The results of NnB{t) are shown in Fig.. [2b- We used a 
least-squares fit of Eq.(14) to the simulation results of the reactive flux, C#s(t) and N HB (t) 
to produce the forward and backward rate constants. We performed the fitting in the short 
time region < t < 0.15 ps to obtain the rate constants and the corresponding average 
hydrogen bond lifetime for the librational relaxation and we also carried out the fitting on 
the longer time region 0.25 ps < t < 3 ps to calculate these quantities for slower rotational 
and translational diffusional relaxation. The inverses of the corresponding forward rate 
constants, which correspond to the average hydrogen bond lifetimes and which we denote 
as l/k short and l/fci ng, are included in Table VII. We note that the value of l//c s hort is very 
close to average hydrogen bond lifetime thb obtained from the continuous hydrogen bond 
correlation function SnB{t) which is not unexpected because both SnB{t) and the short- 
time part of the 'reactive flux' captures the hydrogen bond 'breaking' dynamics due to fast 
librational motion. The hydrogen bond lifetime as given by l/ki ong is significantly longer 
because it excludes the fast librational relaxation. 



VI. SUMMARY AND CONCLUSIONS 



We have investigated the properties of ammonia from the gas phase dimer to the liquid 
state by means of wavefunction based quantum chemistry techniques, density functional 
theory, and ab initio molecular dynamics. The ammonia dimer is a particularly interesting 
hydrogen-bonded system, as it is known to feature a strongly bent hydrogen bond. In order 
to investigate possible consequences of this nonlinear hydrogen bond, important points on 
the potential energy surface were first investigated by means of ab initio methods up to W2- 
extrapolated CCSD(T) theory, see Fig. [T] for the corresponding structures. The structure 
and energetics of important isomers are obtained up to unprecedented accuracy without 
resorting to experiment; see Fig. and Table HVl for structures, and Fig. |2] and Table ITTTI 
for energies. It is confirmed that the structure of the global minimum (having eclipsed C s 
symmetry) possesses a substantially nonlinear hydrogen bond with an H-N- ■ -N angle that 
is predicted to be 20.7°; note that the corresponding angle is predicted to be 5.5° in the 
water dimer. The energy difference between the global minimum and the transition 
state (which is the famous "cyclic structure") is about 3.5 cm -1 , and the energy penalty to 
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make the hydrogen bond artificially linear amounts to only 70 cm -1 , i.e. 0.2 kcal/mol or 
0.009 eV. This implies that the potential energy surface is very flat. 

However, when investigating the ammonia dimer using a wide variety of available 
density functionals (including important representatives of the GGA, meta-GGA and hybrid 
functional families) it was found that none of the functionals checked by us describes the bent 
hydrogen bond in (NH 3 ) 2 satisfactorily. Typically, the hydrogen bond angle is too small by a 
factor of two and the energy of the cyclic C 2 h structure w.r.t. the eclipsed C s global minimum 
is overestimated by a factor of 10 to 30. Since ammonia is an important hydrogen bonded 
liquid, a density functional was developed for studying this subtle system, with particular 
focus on its applicability for condensed phase simulations in the framework of Car-Parrinello 
molecular dynamics. To this end, a novel GGA-type density functional, HCTH/407+, was 
developed, with special emphasis on the nonlinear hydrogen bond and potential energy 
surface of (NH 3 ) 2 . This functional yields the bent hydrogen-bonded equilibrium structure 
of the ammonia dimer, as well as the correct energetics of the low-lying isomers including 
the energy barrier to a linear hydrogen-bond. 

The performance of the new functional in describing the structural and dynamical 
properties of liquid ammonia was investigated by carrying out Car-Parrinello molecular 
dynamics simulations of the liquid phase. In addition, simulations were carried out using 
a GGA-type functional that is widely used for describing associated liquids, BLYP, and 
the results of the two functionals are compared with experiments wherever available. In 
particular, we focused on the atom-atom radial distribution functions, on structure and 
dynamics of hydrogen bonds, and on diffusion as well as orientational relaxation processes. 
Overall, both functionals are found to describe the structural and dynamical properties of the 
liquid phase reasonably well. Importantly, it is shown that the propensity to form a strongly 
bent hydrogen bond - which is characteristic for the equilibrium structure of the gas-phase 
ammonia dimer - is overwhelmed by steric packing effects that clearly dominate the solvation 
shell structure in the liquid state. Similarly, the propensity of ammonia molecules to form 
bifurcated and multifurcated hydrogen bonds in the liquid phase is found to be negligibly 
small. Thus, even functionals that lead to unreasonably linear hydrogen bonds in the 
limiting case of the in vacuo ammonia dimer, such as BLYP, yield a good description of 
liquid ammonia - albeit for the wrong reason! 
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TABLE I: Dissociation energies D c (in kJ/mol) according to Wn the ories for several hydrogen 
bonded dimers compared to the best available estimates from Refs. 



Molecule 


Wl 


W2 


reference 


(HF) 2 


19.13 


19.10 


19.10 ± 0.2 \U7j 


(H 2 0) 2 


20.93 


20.84 


21.1 ± 0.3 [55] 


(H 2 0)(NH 3 ) 


26. .94 


26.82 


27 A ± 3 [118] 
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TABLE II: Relative energies (in cm -1 ) according to Wn theories for the ammonia dimer at the 
respective CCSD(T)/A'PVQZ optimised structures compared to the global minimum; for the latter 
the dissociation energies D e (in kJ/mol) into isolated optimized monomers is reported in the last 
line. In column 2, and column 4, standard Wl and W2 results are reported, while in column 3, a 
counterpoise-correction is employed on each step in the W2 calculation. 



Points on PES 


Wl (w BSSE) 


W2(w/o BSSE) 


W2(w BSSE) 


best estimate 


C s (staggered) 


23.4 




23.5 


23 ± 5 


C2/1 


0.7 


5.9 


2.8 


3.5 ± 3 


linear (eclipsed) 


69.6 




67.5 


67 ± 5 


linear (staggered) 


71.8 




69.9 


70 ± 5 


Global Minimum 


C s (eclipsed) 


13.21 


13.15 


13.13 


13.1 ± 0.2 
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TABLE III: Various first principles results for the dissociation energies D e (in kJ/mol) of the 
eclipsed C s global minimum and relative energies (in cm -1 ) of several structures of the ammonia 
dimer w.r..t. the global minimum. All data, except that for W2, were obtained using an A'PVTZ 
basis set at fully optimised geometries including the counterpoise correction. Wl and our best 
estimates for these values are reported in Table ITT1 



Geometry 
Type 


C s (eel.) 
De 


C s (sta.) 
difference 


C2/1 

difference 


lin (eel.) 
difference 


lin (sta.) 
difference 


Method 


in kJ/mol 


in cm -1 


in cm -1 


in cm -1 


in cm -1 


W2 


13.13 


23.5 


2.8 


67.5 


69.9 


CCSD(T) 


12.22 


18.2 


2.8 


62.8 


63.8 


CCSD 


11.31 


18.4 


3.3 


60.7 


59.4 


MP2 


12.27 


20.3 


2.9 


71.8 


72.7 


HF 


7.58 


12.6 


23.6 


32.3 


34.4 


BLYP 


9.11 


5.7 


97.8 


36.5 


36.4 


PBE 


13.05 


9.2 


91. .5 


42.1 


41.4 


HCTH/120 


9.48 


6.2 


75.1 


31.1 


32.4 


HCTH/407 


11.23 


11.0 


27.0 


39.0 


41.1 


r-HCTH 


10.96 


5.2 


141.8 


29.9 


29.6 


B3LYP 


10.18 


3.3 


76.6 


40.2 


40.1 


B97-1 


12.73 


11.2 


61.1 


41.5 


41.6 


r-HCTH hyb. 


11.27 


9.4 


99.9 


37.7 


37.2 


HCTH/407+ 


13.18 


18.0 


4.0 


51.5 


53.6 
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TABLE IV: Various first principles results for representative structural parameters of the 
dimer, see Fig. ^ using the A'PVTZ basis set. 



Geometry 
Type 


C s (eel.) 

N- • H Distance 


C2/1 

N- • H Distance 


C s (eel.) 
HNN Angle 


C s (sta.) 
HNN Angle 


best estimate 


2.31 


2.52 


20.7 


13.2 


CCSD(T) 1 


2.302 


2.522 


19.86 


13.09 


CCSD(T) 


2.294 


2.527 


16.40 


12.43 


CCSD 


2.331 


2.562 


16.76 


12.39 


MP 2 


2.286 


2. .520 


17.32 


13.28 


HF 


2.541 


2. .772 


13.41 


9.69 


BLYP 


2.341 


2.610 


9.73 


8.55 


PBE 


2.249 


2.512 


9.90 


8.47 


HCTH/120 


2.427 


2.736 


10.27 


8.65 


HCTH/407 


2.493 


2.773 


13.71 


10.25 


r-HCTH 


2.311 


2.645 


8.40 


7.45 


B3LYP 


2.322 


2.573 


10.50 


8.89 


B97-1 


2.298 


2.555 


10.90 


9.00 


r-HCTH hyb. 


2.269 


2.549 


9.57 


8.24 


HCTH/407+ 


2.493 


2.754 


16.95 


11.33 



A'PVQZ basis set 
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TABLE V: Expansion coefficients of the HCTH/407+ functional compared to the coefficients of 
the HCTH/407 functional. 



Functional 


HCTH/407+ 


HCTH/407 


Cl = C X(7 ,0 


1.08018 


1.08184 


C2 = CCaafi 


0.80302 


1.. 18777 


C3 = CCa/3,0 


0.73604 


0.58908 


C4 = C Xa ,l 


-0.4117 


-0..5183 


C5 = CCW,1 


-1.0479 


-2.4029 


C6 = CCa/3,1 


3.0270 


4.4237 


Cl = C X a,2 


2.4368 


3.4256 


C8 = CCW,2 


4.9807 


5.6174 


Cg = CCo/3,2 


-10.075 


-19.222 


ClO = CXa,3 


1.3890 


-2.6290 


Cll = Cc CT( j,3 


-12.890 


-9.1792 


Cl2 = CCa/3,3 


20.611 


42.572 


Cl3 = CXo-,4 


-1.3529 


2.2886 


Cl4 = CCcro-,4 


9.6446 


6.2480 


Cl5 = CCa/3,4 


-29.418 


-42.005 
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TABLE VI: Comparison of the errors of the 407+ fit set for the HCTH/407+ and previous published functional. Also included are the 
errors for properties of nine hydrogen-bonded systems (see text).. 





407+ Set of molecules 


Hydrogen-Bonded Systems 


Property 


RMS Energy 


J2 Gradient 


RMS Dissociation Energy 


RMS H-Bond Shift 


RMS Frequency Shift 


Unit 


[kcal/mol] 


[a.u.] 


[%] 


[%] 


[%} 


HCTH/120 


9.1 


11.51 


8.5 


26.0 


29.1 


HCTH/407 


7.8 


11.91 


7.5 


16.7 


14.8 


HCTH/407+ 


8.0 


11.96 


10.3 


16.8 


12.8 


B3LYP 


9.4 


12.00 


11.0 


30.2 


35.0 


BLYP 


9.7 


19.30 


16.4 


42.1 


42.9 


BP86 


16.4 


17.05 


22.3 


81.9 


68.3 



TABLE VII: Dynamical properties of liquid ammonia. The diffusion coefficient, relaxation times 
and the inverse rate constants are expressed in units of ICT 5 cm 2 sec _1 , ps and ps, respectively. 
Note that all numbers given refer to the fully deuterated species ND3. 



Quantity 


HCTH/407+ 


BLYP 


Experiment 


D 


5.23 


5.44 


8.75 a , 5.95 6 


^_dipole 


0.52 


0.85 




t-NH 


0.37 


0.58 




~HH 
T l 


0.34 


0.54 




dipole 
T 2 


0.18 


0.27 




t-NH 
T 2 


0.15 


0.22 




~HH 

T 2 


0.15 


0.22 




exp 
r 2 






0.26 a , 0.35 6 


TUB 


0.07 


0.10 




1/^short 


0.10 


0.16 




1/Mong 


0.62 


0.60 





at 275 K, b at 252 K 
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Fig. 1. The five relevant optimized structures of the ammonia dimer; dotted lines are only 
guides to the eye.. 

Fig. 2. The potential energy profile obtained from W2 Theory (accurate reference 
calculations) in comparison to that obtained from HCTH/407+ and BLYP functionals at 
the five optimized structures, see insets, from Fig. ^ Note that the profiles obtained from 
other GGA functionals such as PBE are qualitatively similar to that of BLYP. 

Fie. 3. The atom-atom radial distribution functions. The s olid, dashed and the dotted curves 
show the HCTH/407+, BLYP, and the experimental Il2l| results, respectively. 

Fig. 4. The fraction P(n) of molecules having a certain number of accepted n acceptor and 
donated n donor hydrogen bonds, see panels (a) and (b), respectively. Panel (c) shows 
the fraction of hydrogen atoms that donates a certain number of hydrogen bonds. The 
squares and the triangles represent the results for the HCTH/407+ and BLYP functionals, 
respectively, and the dashed lines are drawn to guide the eye. 

Fig. 5. The distribution of the cosine of the N- • -N-H angle for H atoms that belong to the 
nearest neighbours. Panel (a) shows the distribution for hydrogen bonded H atoms (i.e.. 
Ri NH ) < 2.7 A) and panels (b) and (c) show, respectively, the results for non-hydrogen- 
bonded H atoms that appear in the region 2.7 A < R^ NH ^ < 3.7 A and 3.7 A < R^ NH ^ 
< 4.7 A. The solid and the dashed curves show the results for the HCTH/407+ and BLYP 
functionals, respectively. 

Fig. 6. The time dependence various orientational correlation functions, see text for 
definitions and further details. The solid and the dashed curves show the results for the 
HCTH/407+ and BLYP functionals, respectively. 

Fig. 7. The time dependence of various hydrogen bond correlation functions, see text for 
definitions and further details. The different curves are labelled as in Fig.6. 
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